#################################
##' Code for:
##' Title: "Effective climate clubs require ambition, leverage, and insulation"
##' Author: Sam S. Rowan
##' Contact: sam.rowan [at] concordia.ca
##' Date: 2023-10-19
#################################
##' This script:
##' Plotting ideal type clubs
#################################


#### Initialize environment ####

## Packages
library(tidyverse) 
library(countrycode)
library(WDI)

## Working directory
# setwd("/Users/srowa/Dropbox/Projects/climate_clubs_trade")
# -- Set working directory to replication folder

## Set plot theme for later
theme_clubs <-  theme(legend.position = "bottom", 
                      panel.background = element_rect(fill=NA),
                      panel.grid.major.x = element_blank(),
                      panel.grid.major.y = element_blank(),
                      panel.grid.minor.y =  element_blank(),
                      panel.grid.minor.x =  element_blank(),
                      panel.ontop = FALSE,
                      axis.line = element_line(linewidth = 1, colour = "black"),
                      axis.text = element_text(size = 16),
                      legend.text = element_text(size = 12),
                      legend.title = element_text(size = 16),
                      axis.title = element_text(size = 16))


#### Load in covariates for plotting ####

## Not run --- load from disk below
## Load in emissions data from World Development Indicators
# wdi_emissions <- WDI::WDI(indicator = c("SP.POP.TOTL",
#                                         "NE.EXP.GNFS.ZS",
#                                         "EN.ATM.GHGT.KT.CE"),
#                           start = 2014,
#                           end = 2014,
#                           extra = T) %>%
#   filter(income!="Aggregates") %>% 
#   select(iso3c, year,
#          population = SP.POP.TOTL,
#          exports = NE.EXP.GNFS.ZS,
#          ghg_2014 = EN.ATM.GHGT.KT.CE) %>%
#   mutate(ghg_percapita = ghg_2014/population,
#          ghg_multiplicative = ghg_2014 * ghg_percapita,
#          ghg_multiplicative_log = log(1+ghg_multiplicative)) %>%
#   arrange(iso3c) %>%
#   as.data.frame()
# write_csv(wdi_emissions, "data/wdi_emissions.csv")

## Load WDI data from disk
wdi_emissions <- read_csv("data/wdi_emissions.csv")

#### Load in trade data from disk ####
## (made with 'dd_comtrade_bulk_tidy.R') 

comtrade_total <- read_csv("data/comtrade_total.csv")


#### Start by making the clubs ####

## Define EU-27 countries

eu_countries <- c(
  "AUT", # Austria
  "BEL", # Belgium
  "BGR", # Bulgaria
  "HRV", # Croatia
  "CYP", # Cyprus
  "CZE", # Czech Republic
  "DNK", # Denmark
  "EST", # Estonia
  "FIN", # Finland
  "FRA", # France
  "DEU", # Germany
  "GRC", # Greece
  "HUN", # Hungary
  "IRL", # Ireland
  "ITA", # Italy
  "LVA", # Latvia
  "LTU", # Lithuania
  "LUX", # Luxembourg
  "MLT", # Malta
  "NLD", # Netherlands
  "POL", # Poland
  "PRT", # Portugal
  "ROU", # Romania
  "SVK", # Slovakia
  "SVN", # Slovenia
  "ESP", # Spain
  "SWE"  # Sweden
)

## Binary of membership in EU club
eu_club <- wdi_emissions %>% 
  mutate(eu_club = ifelse(iso3c %in% eu_countries, 1, 0)) %>% 
  select(iso3c, eu_club)
# Czechia drops out; no WDI data

## NDC club -- from Robiou du Pont et al. 2017; Rowan 2021
ndc_countries <- c(
  "BGD", # Bangladesh
  "BEN", # Benin
  "BTN", # Bhutan
  "KHM", # Cambodia
  "CMR", # Cameroon
  "TCD", # Chad
  "COM", # Comoros
  "COG", # Congo - Brazzaville
  "COD", # Congo - Kinshasa
  "CIV", # Côte d’Ivoire
  "CYP", # Cyprus
  "DJI", # Djibouti
  "GNQ", # Equatorial Guinea
  "ERI", # Eritrea
  "ETH", # Ethiopia
  "FJI", # Fiji
  "GHA", # Ghana
  "GTM", # Guatemala
  "GIN", # Guinea
  "ISL", # Iceland
  "KEN", # Kenya
  "LBR", # Liberia
  "MDG", # Madagascar
  "MLI", # Mali
  "MLT", # Malta
  "MDA", # Moldova
  "MAR", # Morocco
  "PER", # Peru
  "STP", # São Tomé & Príncipe
  "SWE",  # Sweden
  "CHE", # Switzerland
  "TJK", # Tajikistan
  "TZA", # Tanzania
  "TGO", # Togo
  "UGA"  # Uganda
)


## Binary of membership in NDC club
ndc_club <- wdi_emissions %>% 
  mutate(ndc_club = ifelse(iso3c %in% ndc_countries, 1, 0)) %>% 
  select(iso3c, ndc_club)


## Identifying states with the highest leverage on key emitters

## Get 10 biggest emitters from multiplicative index
biggest_emitters <- wdi_emissions %>% 
  arrange(-ghg_multiplicative) %>% 
  slice(1:10) %>% 
  pull(iso3c)

## Check their names
sort(countrycode(biggest_emitters, 'iso3c', 'country.name'))

## Identify the top 5 markets
hold_export_leverage_10 <- as.vector(NA)
for (emitter in biggest_emitters){
  # print(emitter)
  hold_exporters <- comtrade_total %>% 
    filter(reporter_iso == emitter) %>% 
    arrange(-export_share_to_partner) %>% 
    ungroup() %>% 
    slice(1:5) %>% 
    select(partner_iso) 
  hold_export_leverage_10 <- rbind(hold_export_leverage_10, hold_exporters$partner_iso)
}

## Re-arrange this as a vector of ISO3C
leverage_countries <- hold_export_leverage_10 %>% 
  as.character() %>%
  unique() %>% 
  as.data.frame() %>%
  rename(iso3c = 1) %>%
  filter(iso3c!="NA") %>%
  filter(iso3c!="HKG") %>%
  filter(!(iso3c %in% biggest_emitters)) %>% 
  arrange(iso3c) %>% 
  mutate(leverage_club = 1)

## Check their names
sort(countrycode(leverage_countries$iso3c, 'iso3c', 'country.name'))

## Binary of membership in leverage club
leverage_club <- wdi_emissions %>% 
  mutate(leverage_club = ifelse(iso3c %in% leverage_countries$iso3c, 1, 0)) %>% 
  select(iso3c, leverage_club)

## Put all of these clubs together
clubs_df <- full_join(eu_club, ndc_club, by = "iso3c") %>% 
  full_join(., leverage_club, by = "iso3c") %>% 
  mutate_at(vars(eu_club:leverage_club), ~ifelse(is.na(.), 0, .)) %>% 
  drop_na()


## Make table for SI
cbind.data.frame(`NDC club` = sort(countrycode(clubs_df$iso3c[clubs_df$ndc_club==1], 
                                               'iso3c', 'country.name')),
                 `Leverage club` = c(sort(countrycode(clubs_df$iso3c[clubs_df$leverage_club==1],
                                                      'iso3c', 'country.name')),
                                     rep("", times = length(clubs_df$ndc_club[clubs_df$ndc_club==1])-
                                           length(clubs_df$leverage_club[clubs_df$leverage_club==1])))) %>% 
  knitr::kable(.,
               booktabs = T, 
               linesep = "", 
               format = 'latex')


#### Start with total trade #### 

## Merge all club memberships to dyadic trade data: reporter and partner
total_clubs <- full_join(comtrade_total, clubs_df,
                         by = c("reporter_iso" = "iso3c")) %>%
  rename(reporter_in_eu = eu_club,
         reporter_in_ndc = ndc_club,
         reporter_in_leverage = leverage_club) %>% 
  full_join(., clubs_df, 
            by = c("partner_iso" = "iso3c")) %>%
  rename(partner_in_eu = eu_club,
         partner_in_ndc = ndc_club,
         partner_in_leverage = leverage_club) %>% 
  select(-starts_with("ccode")) %>% 
  mutate_at(vars(reporter_in_eu:partner_in_leverage), ~as.numeric(.)) %>%
  mutate_at(vars(reporter_in_eu:partner_in_leverage), ~ifelse(is.na(.), 0, .))


## Make NDC club and calculate each state's trade exposure
total_ndc_exposure <- total_clubs %>% 
  rename(reporter_in_club = reporter_in_ndc, 
         partner_in_club = partner_in_ndc) %>% 
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1, 
                                   export_share_to_partner, 0)) %>% 
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>% 
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>% 
  ungroup() %>% 
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club) %>%
  distinct() %>% 
  as.data.frame()

# Make leverage club and calculate each state's trade exposure
total_leverage_exposure <- total_clubs %>% 
  rename(reporter_in_club = reporter_in_leverage, 
         partner_in_club = partner_in_leverage) %>% 
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1, 
                                   export_share_to_partner, 0)) %>% 
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>% 
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>% 
  ungroup() %>% 
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club) %>%
  distinct() %>% 
  as.data.frame()

## Make EU club and calculate each state's trade exposure
total_eu_exposure <- total_clubs %>% 
  rename(reporter_in_club = reporter_in_eu, 
         partner_in_club = partner_in_eu) %>% 
  # Create a dummy for whether an export relationship is into the club
  mutate(exporting_toclub = ifelse(partner_in_club==1, 
                                   export_share_to_partner, 0)) %>% 
  # Group to the reporter-club level and sum total trade into the club
  group_by(reporter_iso, partner_in_club) %>% 
  mutate(export_exposure_to_club = sum(export_share_to_partner, na.rm = T)) %>% 
  ungroup() %>% 
  # Retain trade flows into the club
  filter(partner_in_club==1) %>% 
  # Keep only key variables and drop duplicates
  select(reporter_iso, reporter_in_club, export_exposure_to_club) %>%
  distinct() %>% 
  as.data.frame()

## EU leverage over AUS, CHN, USA
# total_eu_exposure %>% filter(reporter_iso %in% c("AUS", "CHN", "USA"))
# total_leverage_exposure %>% filter(reporter_iso %in% biggest_emitters) %>% summarize(median(export_exposure_to_club))


#### Plot membership and trade flows in each club ####

## Merge GHG data into the club trade exposure data

## NDC club: scatterplot
p_ndc1 <- total_ndc_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(x = ghg_multiplicative_log, y = export_exposure_to_club,
                label = reporter_iso, color = reporter_in_club)) +
  geom_smooth(method = "lm", se = F) +
  geom_text(size = 5) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  # scale_y_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_color_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(y = "Export exposure to club members",
       color = "Club member",
       x = expression(paste("GHG emissions" %*% "GHG emissions per capita, log transformed")))

p_ndc2 <- total_ndc_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(export_exposure_to_club,
                group = reporter_in_club,
                fill = reporter_in_club)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(breaks = c(-0.2, 0.2), labels = c("Member", "Non-member")) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  # scale_x_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_fill_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(x = "",
       fill = "Club member",
       y = "Club member")

cowplot::plot_grid(p_ndc1, p_ndc2, align = "v", ncol = 2, rel_widths = c(2/3, 1/3))
ggsave("text/comtrade_ndc_stacked.pdf", units = "in", height = 6, width = 12)

## Leverage club: scatterplot
p_leverage1 <- total_leverage_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(x = ghg_multiplicative_log, y = export_exposure_to_club,
                color = reporter_in_club)) +
  geom_smooth(method = "lm", se = F) +
  geom_text(size = 5, aes(label = reporter_iso)) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  # scale_y_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_color_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(y = "Export exposure to club members",
       color = "Club member",
       x = expression(paste("GHG emissions" %*% "GHG emissions per capita, log transformed"))) 

p_leverage2 <- total_leverage_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(export_exposure_to_club,
                group = reporter_in_club,
                fill = reporter_in_club)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(breaks = c(-0.2, 0.2), labels = c("Member", "Non-member")) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  # scale_x_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_fill_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(x = "",
       fill = "Club member",
       y = "Club member")

cowplot::plot_grid(p_leverage1, p_leverage2, align = "v", ncol = 2, rel_widths = c(2/3, 1/3))
ggsave("text/comtrade_leverage_stacked.pdf", units = "in", height = 6, width = 12)


## For EU club: scatterplot
p_eu1 <- total_eu_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(x = ghg_multiplicative_log, y = export_exposure_to_club,
                label = reporter_iso, color = reporter_in_club)) +
  geom_smooth(method = "lm", se = F) +
  geom_text(size = 5) +
  scale_y_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  # scale_y_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_color_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(y = "Export exposure to club members",
       color = "Club member",
       x = expression(paste("GHG emissions" %*% "GHG emissions per capita, log transformed"))) 

## For EU club: boxplot
p_eu2 <- total_eu_exposure %>% 
  full_join(., wdi_emissions, by = c("reporter_iso" = "iso3c")) %>% 
  mutate(reporter_in_club = as.factor(case_when(reporter_in_club==1 ~ "Member",
                                                reporter_in_club!=1 ~ "Non-member"))) %>% 
  ggplot(., aes(export_exposure_to_club,
                group = reporter_in_club,
                fill = reporter_in_club)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(breaks = c(-0.2, 0.2), labels = c("Member", "Non-member")) +
  scale_x_continuous(limits = c(0,1), breaks = seq(0, 1, 0.1)) +
  # scale_x_continuous(limits = c(0,0.5), breaks = seq(0, .5, 0.1)) +
  scale_fill_manual(values = c("#556B2F", "grey70")) +
  theme_clubs +
  theme(legend.position = "none") +
  labs(x = "",
       fill = "Club member",
       y = "Club member")

cowplot::plot_grid(p_eu1, p_eu2, align = "v", ncol = 2, rel_widths = c(2/3, 1/3))
ggsave("text/comtrade_eu_stacked.pdf", units = "in", height = 6, width = 12)



